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SIGNIFICANCE  AND  EXPLANATION 


Solving  differential  equations  numerically  on  a  computer  is  not  as  easy 
as  it  seems  to  be<  It  is  not  enough  just  to  write  a  computer  program  for  the 
right  hand  side  of  the  differential  equation,  define  the  additional  data  and 
then  call  for  a  solving  subroutine  (e.g.:  Runge  Kutta ,  Adams-Stormer ,  Picard- 
LindelOf).  The  numerical  results  obtained  in  this  way  are  usually  disastrous 
unless  one  selects  the  "correct"  step  size  and  a  "corresponding"  convergence 
order . 

Everybody  working  in  numerical  computation  knows  of  this.  But  it  is 
hardly  ever  mentioned  in  detail  (if  at  all)  in  textbooks. 

What  most  people  do  not  know,  however,  is  the  fact  that  this  problem  can 
be  solved  completely  by  using  interval  methods.  This  is  one  of  the  reasons 
why  this  paper  was  written.  The  above  mentioned  problem  is  demonstrated  and 
resolved  by  using  a  very  simple  example. 

Unfortunately  in  interval  methods  a  new  problem  arises.  It  is  sometimes 
called  the  "wrapping  effect"  and  can  produce  an  unwanted  explosion  of  the 
computed  error  bounds.  This  effect  is  also  explained  in  this  paper  and  a 
solution  to  it  is  presented. 

Finally  a  survey  is  given  to  show  all  the  possible  numerical  interval 
methods  which  are  known  up  to  now.  It  is  pointed  out  that  as  a  byproduct  of 
using  interval  methods  the  problem  of  stability  just  vanishes,  if  one  uses  the 
appropiate  interval  termination  criterium.  The  list  of  references  given  here 
is,  with  123  entries,  very  long  and  includes  all  the  relevant  work  in  this 
field  which  is  known  to  the  author. 


USING  INTERVAL  METHODS  FOR  THE  NUMERICAL  SOLUTION  OF  ODE'S  ) 


Karl  L.E.  Nickel 

1 . Introduction 

In  this  paper  the  initial  value  problem 


(1) 

u'(t)  -  f(t,u(t)) 

for  t  e  I  , 

(2) 

u{  a)  -  a 

is  treated.  Here  a,b  are  real  numbers  with  a  <  b,  and  I  im  [a,b] .  Let  n 
be  an  integer  and  ae  tf1,  u  s  I  +  rf1,  £  :  I  x  rf1  tf1.  The  functions  f 
and  u  have  to  satisfy  certain  regularity  conditions. 

A 

Let  u(t)  be  a  solution  of  (1),  (2).  Such  a  solution  is  not  explicitly 
known  in  general.  Also,  it  does  not  have  to  be  uniquely  determined.  Wanted 
are  explicitly  computable  error  bounding  functions 
v,w  t  I  ♦  rf1 

A 

such  that  v  and  w  are  respectively  lower  and  upper  bounds  for  u,  i.e. 

A 

(3)  ®v(t)  <  u(t)  <  w(t)  for  tel. 

Here  and  in  what  follows,  the  order  relation  <  is  defined  componentwise.  I 
the  space  of  the  functions  v,w  the  symbol  <  means  the  correspondingly 
induced  partial  ordering.  Using  the  notation  [v,w]  for  the  order  interval 
of  the  factions  v,w,  the  inequality  (3)  can  be  regarded  as  the  inclusion 


1 )  This  paper  is  a  condensed  version  of  a  main  lecture  given  at  the 
Second  Seminar  On  The  Numerical  Treatment  Of  Differential  Equations 
in  Halle/Saale,  DDR  (East  Germany)  at  May  16  to  19  of  1983. 


(4) 


u  e  [v,w] 


Suppose  the  problem  of  finding  bounds  v,w  satisfying  (3)  or  (4)  has  been 

solved.  Then  It  Is  quite  natural  to  ask  for  bounds  v,w  to  all  solutions 
* 

u(t)  of  the  same  differential  equation  (1),  but  with  different  Initial 
values.  Let  A  _c_  sP  be  a  bounded  set.  Then  one  can  replace  the  initial 
condition  (2)  by  the  initial  inclusion 

(5)  u(a)  e  A 

and  ask  for  bounds  (3)  or  (4)  to  all  solutions  of  (1),  (5)  . 

In  general  it  is  quite  easy  to  write  down  rough  bounds  which  satisfy  (3) 
or  (4).  It  is  not  as  simple,  however,  to  find  bounds  which  are  "small”, 
"realistic"  or  even  "optimal".  Define  the  "interval  hull"  of  the  set 

a  a 

{u}  of  all  solutions  u  of  (1)  by 

A  A  A 

(6)  hull  {u}  s«  [inf  {u},  sup  {u}]  . 

A  A 

If  (1),  (2)  has  a  unique  solution  then  obviously  hull  {u}  *  u.  The  beyt  one 

A 

may  hope  to  achieve  is  to  determine  v,w  such  that  [v,w]  *  hull  {u} .  There 

A 

are  theoretical  investigations  of  the  evaluation  of  hull  {u}  ,  see  Nickel 
[79]  and  [80]  .  Due  to  the  limitations  of  numerical  computations,  one 
normally  has  to  settle  for  an  (outer)  approximation  of 

A 

hull  (u)  • 


2.  Why  Computing  Interval  Bounds  Instead  of  Approximations  ? 


It  has  been  exactly  25  years  ago  since  the  first  electronic  computer  was 
installed  at  the  University  of  Karlsruhe/GERMANY  in  the  summer  semester  of 
1958.  The  author  of  this  paper  was  given  the  responsibility  for  this  machine 
(  Zuse  Z  22).  One  of  the  very  first  problems  I  treated  numerically  with  it 
was  the  solution  of  the  following  initial  value  problem  (n  «  1): 

(7)  u'  -  |1-u2|,  u(0)  -  0  . 

The  uniquely  determined  solution  of  (7)  is  obviously  u(t)  ;*  tanh  t.  It  has 

4k 

the  property  0  <  u  <  1.  Some  of  the  results  obtained  with  a  Runge-Kutta 
method  are  sketched  in  Figure  1.  The  four  different  step  sizes  h  :=»  0.85, 

0.90,  0.95,  1.0  were  used.  The  values  actually  computed  are  denoted  by 
circles.  The  lines  between  them  are  dram  to  link  up  the  results  with  the 
same  step  size. 


Figure  1 .  Numerical  solution  of  the  initial  value  problem  (7)  by  using  a 

Runge-Kutta  method  with  the  four  step  sizes  h  :■  0.85,  0.90,  0.95,  1.0. 


The  numerical  result  displayed  in  Figure  1  is  rather  embarassing  and 
displeasing:  three  out  of  four  examples  enter  the  "forbidden”  region  u  >  1 
and  there  grow  exponentially  fast  toward  +  *  .  Furthermore,  the  dependence 
upon  the  step  size  parameter  h  is  not  monotone,  as  one  would  expect. 

If  one  solves  this  problem  (7)  using  many  different  methods  with  many 
different  step  sizes,  the  following  occurs  (see  Figure  2)  :  One  sees  a  cluster 

A 

of  circles  without  any  indication  of  the  whereabouts  of  the  solution  u.  Note 
that  the  values  high  above  the  "critical  line"  u  •  1  belong  to  both  very 
large  and  very  small  step  sizes  h  (1)  .  As  bad  as  this  is,  worse  is  to 
come:  By  adding  the  numerical  results  of  one  more  method  or  one  more  step 

A 

size  no  additional  information  about  the  solution  u  can  be  gained  ! 


Figure  2.  Sketch  of  numerical  results  of  solving  the  Initial  value 
problem  (7)  by  using  different  methods  and  different  step  sizes 
This  sketch  is  schematic,  not  actually  computed. 


Por  the  "normal"  user  of  a  computer  and  of  the  numerical  methods  praised 


in  mathematical  textbooks  such  a  situation  is  a  nightmare  and  most  certainly 

intolerable.  By  doing  calculations  over  and  over  again,  nothing  is  learned 

about  the  values  of  the  desired  solution  u(t).  I  myself  have,  therefore,  used 

this  special  example  (7)  and  Figure  1  in  classes  over  the  years  again  and 

again  to  teach  my  students  a  sound  distrust  of  the  "usual"  textbook  methods. 

By  using  interval  methods  for  solving  (7)  the  whole  conception  changes. 

This  is  illustrated  in  Figure  3.  An  interval  method  produces  bound  functions 

v,w  such  that  (4)  is  valid.  A  first  very  rough  inclusion  to  the  solution  of 
*  -t 

(7)  is  u(t)  e  [v,w]  :=  [1  -  e  ,  1]  .  This  can  be  found  easily  with  the  aid 
of  differential  inequalities  (see  W.  Walter  [117]).  A  second  inclusion  is 
computed  by  using  a  power  series  method  (see  Section  5)  on  a  coarse  grid  with 
the  step  size  h  :*  0.5.  These  two  interval  bounds  are  combined  in  Figure  3. 
The  first  one  is  realistic  for  large  values  of  t,  while  the  second  gives 


very  pessimistic  results.  For  0  _<  t  <.  1  the  opposite  is  true.  There  are 
two  very  pleasant  and  surprising  results  : 

i)  The  bounds  are  valid  for  all  values  of  t,  not  only  at  the  grid 
points. 

ii)  The  intersection  of  the  two  inaccurate  interval  inclusions  gives 
a  much  better  result  in  the  whole  line  0  <  t  <  •  than  each  one 
produces  alone. 

Hence,  by  taking  into  account  the  results  of  more  and  more  interval 

A 

computations,  one  gets  more  and  more  information  about  the  solution  u.  Th: 
result  is  completely  contrary  to  the  result  gained  in  the  non-interval 
numerical  computation  (see  Figures  1  and  2). 


3.  Blowing-Up  Of  Intervals;  A  Difficulty  Inherent  In  Ml  Computations 


Involving  Seta  . 

There  is  a  certain  built-in  difficulty  in  solving  initial  value  problems 
with  inclusions  using  general  sets.  Sometimes  it  is  called  the  "wrapping 
effect".  It  has  first  been  noticed  with  interval  methods  by  R.  B.  Moore  [69] 
but  this  difficulty  arises  also  even  more  drastically  when  using  norm 
bounds.  It  is  best  explained  by  Moore's  own  example.  Considered  is  the 
system  of  two  linear  differential  systems  (n  -  2) 


The  solutions  are  concentric  circles  around  the  origin  in  the  u^u^  plane. 
Suppose  the  initial  values  for  t  s  0  lie  in  a  (small)  square  A  =  [a]  ,  see 
Figure  4.  Then  for  t  >  0  this  square  is  rotated  around  the  origin. 
Bounding  this  square  by  an  interval  according  to  (3)  and  (4)  has  the 


following  meaning  in  the  ,u2-*plane:  One  has  to  find  a  new  enclosing 

rectangle  (a  square)  with  sides  parallel  to  the  coordinate  axes.  In  Figure  4 

this  is  first  done  at  t  =  x/6,  producing  a  larger  square.  Now  this  new 

and  larger  square  also  rotates  for  t  >  n/6  and  has  to  be  bounded  again  by  a 

rectangle  (square)  parallel  to  the  axes.  By  repeating  this,  the  enclosing 

squares  become  larger  and  larger,  independently  of  the  fact  that  the  image  of 

the  initial  set  [a]  keeps  its  size.  It  can  be  seen  easily  that  the  original 

2w 

square  is  blown  up  by  a  factor  of  e  *  535  (t)  after  only  one  revolution 
(t  *  2ir)  as  the  step  size  h  vanishes.  This  is  the  price  to  be  paid  for 
bounding  the  solutions  of  (8)  with  initial  values  in  [a]  by  rectangles 
parallel  to  the  axes,  i.e.  intervals. 

Obviously  this  problem  with  the  special  system  (8)  could  be  resolved  by 
using  a  disc  to  enclose  the  initial  data  instead  of  the  square  [a]  .  Such  a 
disc  keeps  its  shape  and  size  while  rotating  for  t  >  0.  Generalizing  this 
idea  to  arbitrary  systems  would  mean  replacing  interval  bounds  by  norm  bounds 
(spheres).  Unfortunately,  this  is  even  more  likely  to  produce  exploding 
bounds,  as  can  be  seen  by  very  simple  examples,  even  with  linear  differential 
systems  ( 1 ) . 

In  the  past  15  years  several  ideas  have  been  invented  to  overcome  that 
difficulty: 

i)  Moore  [69]  uses  local  coordinate  transformations. 

ii)  Kahan  [44]  bounds  the  solutions  by  ellipsoids  with  axes  not 
nessessarily  parallel  to  the  coordinate  axes. 

iii)  Walzel  [118]  uses  a  pretransformation  to  get  (hopefully)  a  "better" 
differential  equation. 

iv)  Eijgenraam  [26]  uses  a  solution  set  which  has  built-in  coordinate 


transformations 


Some  years  *<30  R.  Lohner  [59]  and  the  author  [79]  solved  this  problem 
independently  of  each  other  at  least  for  linear  systems  (1),  by  using  the 
following  facts:  Zf  the  right  hand  side  f  in  (1)  is  linear  then  the  mapping 

A  n  *  n 

from  the  initial  value  u(a)  *  a  e  R  to  u(t)  e  R  for  a  fixed  value  of 
t  >  a  is  obviously  an  affine  transformation.  Now  take  any  m  >  n  points 

A 

in  the  u-space  and  let  u^(t)  be  the  solutions  of  the  differential 

A 

equation  (1)  for  the  initial  conditions  u^fa)  ■  0^  for  i  -  1(1 )m. 

Let  A  be  the  set  spanned  by  the  {a^},  i.e.  the  convex  hull  of  the  points 

A 

cx^  for  i  =*  1(1)  m  .  Let  U(t)  be  the  set  spanned  similarily  by  the  {u^(t)} 
for  a  fixed  value  of  t  >  a.  Then,  obviously, 

A 

u(t)  e  U(t)  for  t  e  I 


for  any  solution  u(t)  of  (1)  and  (5),  see  the  Figures  5  and  6.  If  the 


t  -a 


t>a 


Figure  5.  Affine  transformation  of  the  initial  interval  set  A  *  [a] 
to  U(t)  for  t  >  a. 


initial  set  A  is  an  n-dimensional  interval  defined  by  the  2n  "upper"  and 
"lower"  corners,  then  the  transformed  set  U(t)  is  in  general  not  an  interval 
(see  Figure  5).  But  it  can  be  described  completely  by  using  only  the  2n 

A 

points  of  the  values  u^(t)  .  If  one  starts  with  a  simplex  A,  then  U(t)  is 
also  a  simplex  defined  by  ra  =  n  +  1  points,  see  Figure  6.  Hence  for  linear 
differential  equations  it  suffices  to  solve  n+1  initial  value  problems  (1), 

A 

(2)  in  order  to  bound  all  the  infinitely  many  solutions  u(t)  of  (1)  with 


initial  values  (5).  After  getting  such  simplex  bounds ,  it  poses  no  problem  to 


transform  them  into  interval  bounds  by  talcing  hull  U(t)  for  any  tel. 


t  *  <*  t  >  oi 


Figure  6.  Affine  transformation  of  the  initial  simplex  set  A  to  the  new 
simplex  U(t)  for  t  >  a  together  with  the  two  intervals  hull  A  and 
hull  U(t)  (dashed  lines). 


In  general,  the  functions  u^(t)  cannot  be  evaluated  explicitly.  Hence 
they  have  to  be  computed  numerically  by  using  interval  methods.  By  moving 

A 

from  t  to  t+h,  therefore,  a  (small)  interval  containing  u^(t+h)  is 
produced.  Hence  the  "corners"  of  U(t)  have  to  be  redefined  at  each  h-step, 
see  Figure  7.  This  gives  a  "computable"  set  U(t)  instead  of  the  optimal 
set  U(t)  for  which  inclusion  U(t)  c  U(t)  is  satisfied. 


These  ideas  have  been  outlined  and  elaborated  in  a  paper  of  J.  Conradt 
[19].  The  method  worked  well,  beyond  all  expectations.  Results  are  given  in 
the  next  section. 


Figure  7.  The  computable  set  U(t)  (dashed  lines)  containing  U(t). 


4.  A  Posteriori  Bounds 


In  this  Section  it  is  assumed  that  an  approximation  u  to  the  solution 

A 

u  of  (1),  (2)  is  known.  Two  error  bounding  functions  o,  o  :  I  ♦  if*  which 

IV 

are  computed  from  u  are  wanted  such  that  the  following  inclusion  is  true  for 
A  ~ 

the  error  u  -  u  : 

(9)  *  u(t)  -  u(t)  <  a(t)  for  tel. 

To  get  such  a  function  u  :  I  ♦  sP  (which  has  to  be  defined  for  all  values 
of  tel)  one  usually  takes  a  pointwise  approximation  for  a  certain  step  size 

A  IV 

h  and  transforms  it  to  u(t)  by  using  spline  interpolation.  Once  that  u 
is  known,  the  error  bounds  a  ,  a  can  be  computed  by  using  known  theorems 
from  the  theory  of  differential  inequalities  (see  W.  Walter  [117]).  This 
method  has  been  worked  out  first  by  Markowitz  [61],  [62].  Later  Conradt  [19] 
made  improvements,  gave  complete  interval  programs  and  computed  several 
examples . 

A  very  important  property  of  (9)  is  that  it  gives  an  improvement  (!)  of 

a,  • 

the  approximation  u  if  the  two  bounds  0  and  o  both  have  the  same  sign. 

Experience  shows  that  this  can  often  be  accomplished  I 

In  his  paper  Conradt  [19]  treated  the  differential  equation  (8)  as  an 
example  for  his  method.  The  initial  values  u^tO)  s*  0  ,  u2(0)  1  were 

A  A 

chosen.  This  gives  the  solution  u^(t)  *  sin  t  ,  u^(t)  *  cos  t  •  An 

approximation  u  was  evaluated  by  using  a  4th  order  Runge-Kutta  method  and  by 
interpolating  with  4th  order  splines.  The  following  two  results  are  typical. 

They  where  obtained  by  using  the  step  size  h  :■  0.125  with  single  precision 

I 


(■  7  decimal  digits)  on  a  UNIVAC  1110  computer 


i)  Naive  Interval  method  s  At  t  ■  52  he  got  very  pessimistic  bounds 


with 


8-4io+16  * 


li)  Point  enclosure  method  from  Section  3  s  At  t  *  2152  (an  argument, 
more  than  41  times  larger  than  in  i) )  he  found  the  much  better 
values 


2i  * 

1.03i0-3  , 

°1  “ 

2,31 10~3 

IQ 

M 

1 

-4*66io*3  ' 

®2  " 

-3-46io“3 

Note  that  these  two  pairs  of  bounds  0^  ,  and  o2  ,  each  have 
the  same  sign.  Hence  his  computed  bounds  give  a  considerable 
inproveunt  of  hi.  approximation.  S,  ^  111 

In  Figure  8  the  two  strips  o^,  and  are  sketched  for  the 

-  ~  * 

case  ii).  They  contain  the  two  errors  u^  -  u^  and  u2  -  u2  for 
0  <_  t  £  94.  Note  that  the  error  terms  are  periodic  and  change  sign  as  the 
solution  does.  This  Figure  was  taken  from  Conradt  [19]. 


5.  Which  Interval  Methods  Are  Available  To  The  Uaer  ? 

There  are  3  claseee  of  methods  for  proving  the  existence  of  a  solution  to 
the  Initial  value  problem  (1),  (2).  With  any  one  of  these,  a  class  of 
corresponding  numerical  approximation  methods  can  be  associated.  This 
correspondence  is  shorn  in  Table  1. 

Table  1 . 


Existence  Theorem 

Numerical  Method 

Problems  /  Advantages 

Peano 

Finite  Differences 

Stability  ?  /  Any  (high) 

order  of  stability 

possible 

Picard-Lindelflf 

Numerical  Quadrature 

and  Iteration 

Only  linearly  convergent  / 

Always  stable 

Power  Series  Method 

Taylor  Series 

f  has  to  be  analytic  / 

Automatic  generation  of 

coefficients 

In  what  follows  the  expression  "point-method"  will  be  used  to  distinguish 
between  the  usual  methods  producing  points  and  interval  methods  which  produce 
inclusions  (3)  or  (4)  as  results.  During  the  last  years  some  progress  has 
been  made  in  transforming  some  of  the  existing  point-methods  into  interval 
methods  or  in  creating  new  interval  procedures.  This  progress  is,  however. 


still  much  too  slow.  One  reason  for  this  may  be  that  far  too  few  people  are 


working  in  this  field.  If  only  10%  of  the  many  mathematicians  working  with 
numerical  methods  for  solving  ODE  's  would  choose  interval  methods  much  more 
could  be  achieved  in  a  very  few  years. 

5.1  Finite  Differences 

Nearly  no  interval  difference  methods  have  been  published  up  to  now.  One 
of  the  few  exceptions  is  the  paper  of  Shokin  [103]  who  gives  an  interval 
Runge-Kutta  method;  see  also  Oelschlaegel  and  Nitsche  [81]  on  this  paper. 

One  of  the  easiest-to-generate  methods  is  an  interval  extension  of 
Euler's  method  which  is  repeatedly  mentioned  in  the  literature.  Because  of  a 
poor  convergence  rate  (only  linear)  it  is,  however,  hardly  ever  used  for 
practical  purposes. 


5.1,1  Stability 


In  difference  methods  stability  is  the  most  important  problem.  Dahlquist 
and  Rutishauser  detected  this  independently  of  each  other  in  1951.  If  a 
point-method  is  unstable,  then  no  convergence  of  the  approximation  to  the 
solution  can  be  expected.  Surprisingly  enough  the  same  is  not  true  for 
interval  methods  I  This  has  been  shown  by  the  author  [77] ,  provided  that  the 
natural  termination  criterion  for  interval  sequences  is  used.  Hence  with 
interval  methods  one  can  use  unstable  difference  schemes  and  still  have 
convergence  J  The  only  effect  is  that  the  rate  of  convergence  will  slow 
down.  More  precisely  :  In  interval  computation  combined  with  a  termination 
criterion  the  order  of  convergence  is  the  difference  between  the  order  of 
consistency  and  the  order  of  (in-)  stability. 


15 


5.2  Picard-LindelSf  Iteration 


This  is  one  of  the  "classic"  methods  which  use  Interval  Analysis.  The 
two  equations  (1)/  (2)  are  written  as  one  Vol terra  integral  equation  and  then 
solved  iteratively  by 

*  uQ(t)  arbitrary  , 

J 

'  Uv+1(t)  :*  o  +  /  f(s,uv(s))  ds  for  v  »  0,  1,  ...  . 

0 

In  interval  methods  all  the  functions  uy,  f  in  (10)  are  regarded  as  interval 
functions  (interval  extensions  of  the  point-functions).  There  are  three  steps 
to  be  taken  : 

A 

i)  Find  a  priori  bounds  v,w  for  the  solution  u  of  (1),  (2),  such 

A 

that  u  e  uQ  [v,  w]  .  Sometimes  it  is  not  easy  to  find  close 
bounds.  If  Uq  is  pessimistic,  more  iterations  are  necessary. 

ii)  Find  an  interval  extension  to  f  which  is  "close".  If  f  is 
rational  this  can  be  done  automatically. 

iii)  Perform  the  integration  in  (10)  as  interval  integration  using  one 
of  many  known  procedures. 


One  big  advantage  in  using  (10)  is  that  no  stability  problems  occur.  One 


of  the  disadvantages,  however,  is  the  slow  (linear)  convergence.  This  may  be 


the  reason  why,  up  to  now,  this  class  of  methods  have  been  used  only 
occasionally . 

Quite  recently,  however,  a  new  approach  to  (10)  has  been  published  by 
Raith  184].  Re  generates  a  method  of  arbitrarily  high  (I)  order  of 
convergence.  A  big  advantage  of  his  method  is  that  the  regularity  assumptions 
imposed  on  the  right  hand  side  f  of  (1)  are  very  weak  (continuity  normally 
suffices).  Kxperiments  show  good  results. 

5.3  Power  Series  Method 


Let  f  be  analytic  in  lx  wP.  Ulan  (1),  (2)  has  a  unique  solution  I 

which  is  also  analytic.  Assise  that  it  exists  for  all  tel.  Let  m  >  1  be 
an  integer  number  and  0  <  h  ±  b-a  a  fixed  step  sice.  Then  by  Taylor's 

A 

Theorem,  u(t)  satisfies  the  recurrence  relation 
»  «-1  *  (y )  "(m) 

(11)  u(t+h)  -  l  —  hW  ♦  u  h*  with  t  <  C  <  t+h  . 

y-0  w»  ■* 

Assume  that  intervals  U^(t)  can  be  computed  for  all  tel, 

V  "  0(1)m,  such  that  the  following  inclusions  are  trust 

I 

(12)  u  (M)(t)  e  o  <t)  for  u  -  0(1)  m-1 
and 

(13)  u  (**(C>  e  U  (t)  I 

m 

for  all  t  e  [a,b-h],  t  <  C  <  t+h,  0  <  h  <  b-a. 

I 
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Define  for  simplicity  U(t)  s“  Ug(t)  and  the  grid  t  :■  a  +  vh 
for  v  -  0(  1  )k  with  k-1  <  (b-a)/h  <  k.  Then  by  (11),  (12)  and  (13)  the 
following  interval  method  is  suggested  for  solving  (1),  (5)  recursively: 


for  0  <  t  <  h 


Note  that  U(t)  is  defined  by  (14)  for  all  t  e  1  .  By  (11),  (12),  (13)  and 

(14)  the  following  (desired)  inclusion  is  guaranteed: 

* 

u(t)  e  U(t)  for  all  tel. 

If  all  intervals  U^(ty)  are  globally  bounded,  then  method  (14)  obviously 
defines  a  convergent  sequence  with 

A 

lim  U(t)  “  u(t)  for  all  t  £  I  . 
h-*-0 

It  seems  at  first  that  the  conditions  imposed  on  f  and  u  are  very 
strict  and  that,  moreover,  the  assumptions  concerning  the  computability  of  the 
intervals  by  (12)  and  (13)  can  rarely  be  met.  Unfortunately,  this  wrong 

opinion  is  still  believed  by  most  numerical  analysts.  It  has  been  shown, 
however,  by  people  working  in  interval  analysis  that  the  opposite  is  true  and 


that  the  method  presented  in  (14)  has  nearly  no  limitations  on  it.  There  are 
two  reasons  for  this: 

i)  In  numerical  analysis  any  function  f  has  to  be  evaluated  by  a 
program  on  a  computer.  Any  program,  however,  can  only  produce  the 
values  of  a  piecewise  rational  function,  and  so  any  "usable" 
function  f  is  piecewise  analytic. 

ii)  It  can  be  shown  that  the  computation  of  the  intervals  ( t)  in 

(12)  and  (13)  can  be  carried  out  automatically  by  software  using  the 
compiler  of  the  computer,  if  only  f  is  given  as  a  program  with  no 
juaps.  This  is  known  as  "  automatic  differentiation  ",  see 

R.  E.  Moore  [68],  I>.  B.  Rail  [87]  and  others. 

As  an  Example  the  initial  value  problem  (7)  was  treated  with  that 
method,  using  m  -  5  (convergence  order  m-1  -  4)  and  using  different  step 
sizes  h  i“  1,  1/2,  1/4,  1/8,  1/16,  1/32.  The  results  are  shown  in  Figure  9. 
There  only  the  bounds  are  sketched.  In  order  not  to  overload  the  picture  the 
computed  interval  strips  have  not  been  hatched  as  in  Figure  3.  The  bounds 
shown  in  Figure  3  have  been  computed  with  this  method,  too,  using  the  step 
size  h:«  0.5  . 

This  method  and  Raith's  method  described  in  Section  5.2  both  work 
remarkably  well.  Op  to  now,  no  data  are  available  which  favor  one  method  over 


the  other 


Figure  9.  Solution  of  problem  (7)  by  using  the  power  series  method  with 
the  different  step  sizes  h  »■  1,  1/2,  1/4,  1/8,  1/16,  1/32. 
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